NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY,  CALIFORNIA 


THESIS 


SIMULATING  TSUNAMIS  IN  THE  INDIAN  OCEAN 
WITH  REAL  BATHYMETRY  BY  USING  A  HIGH-ORDER 
TRIANGULAR  DISCONTINUOUS  GALERKIN 
OCEANIC  SHALLOW  WATER  MODEL 

by 

Dimitrios  Alevras 
March  2009 

Thesis  Co-Advisors:  Francis  X.  Giraldo 

Timour  Radko 


Approved  for  public  release;  distribution  is  unlimited 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


j  REPORT  DOCUMENTATION  PAGE 

Form  Approved  OMB  No.  0704-0188  f 

Public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instruction, 
searching  existing  data  sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send 
comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information,  including  suggestions  for  reducing  this  burden,  to 
Washington  headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington,  VA 
22202-4302,  and  to  the  Office  of  Management  and  Budget,  Paperwork  Reduction  Project  (0704-0188)  Washington  DC  20503. 

1.  AGENCY  USE  ONLY  (Leave  blank)  2.  REPORT  DATE  3.  REPORT  TYPE  AND  DATES  COVERED 

March  2009  Master’s  Thesis 

4.  TITLE  AND  SUBTITLE 

Simulating  Tsunamis  in  the  Indian  Ocean  with  Real  Bathymetry  by  using  a  High- 
Order  Triangular  Discontinuous  Galerkin  Oceanic  Shallow  Water  Model. 

5.  FUNDING  NUMBERS 

6.  AUTHOR(S) 

Lt  Dimitrios  Alevras,  Hellenic  Navy 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Naval  Postgraduate  School 

Monterey,  CA  93943-5000 

8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

9.  SPONSORING  /MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

N/A 

10.  SPONSORING/MONITORING 
AGENCY  REPORT  NUMBER 

11.  SUPPLEMENTARY  NOTES  The  views  expressed  in  this  thesis  are  those  of  the  author  and  do  not  reflect  the  official  policy 
or  position  of  the  Department  of  Defense  or  the  U.S.  Government. 

12a.  DISTRIBUTION  /  AVAILABILITY  STATEMENT 

Approved  for  public  release;  distribution  is  unlimited. 

12b.  DISTRIBUTION  CODE 

13.  ABSTRACT  (maximum  200  words) 

The  discontinuous  Galerkin  (DG)  method  has  been  accepted  in  the  last  decade  by  the  geosciences  community  as  an 
important  component  of  geophysical  fluid  dynamics.  The  high-order  accuracy,  geometric  flexibility  to  use 
unstructured  grids,  local  conservation,  and  monotonicity  properties  of  the  DG  method  make  it  a  prime  candidate  for 
the  construction  of  future  ocean  and  shallow  water  models. 

This  study  focuses  on  formatting  real  bathymetry  data  of  the  Indian  Ocean  in  order  to  simulate  the  propagation 
stage  of  the  Indian  Ocean  tsunami  that  occurred  on  December  26,  2004,  by  using  a  DG  model.  In  order  to  validate  this 
simulation  the  study  uses  real  measurements.  The  model  results  are  compared  to  tide  gauge  data  from  several  stations 
around  the  Indian  Ocean,  satellite  altimetry,  and  field  measurements.  These  results  show  that  the  model  gives  accurate 
estimates  of  arrival  times  in  distant  locations. 

14.  SUBJECT  TERMS  Tsunami  Simulation,  Indian  Ocean  Tsunami  2004,  Triangular 
Discontinuous  Galerkin  Method,  Propagation  stage,  Oceanic  Shallow  Water  Model 

15.  NUMBER  OF 

PAGES 

113 

16.  PRICE  CODE 

17.  SECURITY 
CLASSIFICATION  OF 
REPORT 

Unclassified 

18.  SECURITY 

CLASSIFICATION  OF  THIS 
PAGE 

Unclassified 

19.  SECURITY 
CLASSIFICATION  OF 
ABSTRACT 

Unclassified 

20.  LIMITATION  OF 
ABSTRACT 

UU 

NSN  7540-0 1  -280-5500  Standard  Form  298  (Rev.  2-89) 


Prescribed  by  ANSI  Std.  239-18 


1 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


11 


Approved  for  public  release;  distribution  is  unlimited 


SIMULATING  TSUNAMIS  IN  THE  INDIAN  OCEAN  WITH  REAL 
BATHYMETRY  BY  USING  A  HIGH-ORDER  TRIANGULAR 
DISCONTINUOUS  GALERKIN  OCEANIC  SHALLOW  WATER  MODEL 


Dimitrios  Alevras 
Lieutenant,  Hellenic  Navy 
B.A.,  Hellenic  Naval  Academy,  1994 


Submitted  in  partial  fulfillment  of  the 
requirements  for  the  degree  of 


MASTER  OF  SCIENCE  IN  METEOROLOGY  AND  PHYSICAL 

OCEANOGRAPHY 


And 

MASTER  OF  SCIENCE  IN  APPLIED  MATHEMATICS 


from  the 


NAVAL  POSTGRADUATE  SCHOOL 
March  2009 


Author:  Dimitrios  Alevras 


Approved  by:  Francis  X.  Giraldo, 

Co-Advisor 

Timour  Radko, 

Co-Advisor 

Carlos  F.  Borges,  Chainnan 
Department  of  Applied  Mathematics 

Jeffrey  D.  Paduan,  Chainnan 
Department  of  Oceanography 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


IV 


ABSTRACT 


The  discontinuous  Galerkin  (DG)  method  has  been  accepted  in  the  last  decade 
by  the  geosciences  community  as  an  important  component  of  geophysical  fluid  dy¬ 
namics.  The  high-order  accuracy,  geometric  flexibility  to  use  unstructured  grids, 
local  conservation,  and  monotonicity  properties  of  the  DG  method  make  it  a  prime 
candidate  for  the  construction  of  future  ocean  and  shallow  water  models. 

This  study  focuses  on  formatting  real  bathymetry  data  of  the  Indian  Ocean  in 
order  to  simulate  the  propagation  stage  of  the  Indian  Ocean  tsunami  that  occurred 
on  December  26,  2004,  by  using  a  DG  model.  In  order  to  validate  this  simulation 
the  study  uses  real  measurements.  The  model  results  are  compared  to  tide  gauge 
data  from  several  stations  around  the  Indian  Ocean,  satellite  altimetry,  and  held 
measurements.  These  results  show  that  the  model  gives  accurate  estimates  of  arrival 
times  in  distant  locations. 


v 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


vi 


TABLE  OF  CONTENTS 


I.  INTRODUCTION .  1 

II.  BACKGROUND .  5 

A.  GOVERNING  EQUATIONS .  5 

1.  Shallow  Water  Equations .  5 

B.  TRIANGULAR  DISCONTINUOUS  GALERKIN  METHOD  .  .  6 

1.  Basis  Functions .  6 

2.  Integration .  8 

3.  Semi- discrete  Equations  .  8 

4.  Matrix  Form  of  the  Semi-discrete  Equations .  9 

5.  Time  Integrator .  11 

III.  SIMULATION  OF  THE  INDIAN  OCEAN  TSUNAMI .  13 

A.  MODEL  DESCRIPTION  .  13 

1.  Boundary  Conditions .  13 

2.  Tides  /  Coriolis  /  Viscosity .  13 

B.  REAL  BATHYMETRY  DATA  .  14 

C.  TSUNAMI  SOURCE  MODEL .  15 

D.  REAL  MEASUREMENT  DATA  .  17 

1.  Tide  Gauge  Records  .  17 

2.  Satellite  Altimetry  .  25 

3.  Field  Measurements .  27 

IV.  SIMULATIONS  RESULTS .  29 

A.  COMPARISON  TO  TIDE  GAUGE  RECORDS .  29 

1.  Tide  Gauge  Stations  Located  at  the  Tropical  Indian  Ocean 

and  West  of  India  Regions .  29 

2.  Tide  Gauge  Stations  Located  South-southeast  of  India  .  51 

B.  COMPARISON  TO  SATELLITE  ALTIMETRY .  72 

vii 


1.  Jason  1 .  72 

2.  Topex/Poseidon .  74 

3.  Envisat .  76 

C.  COMPARISON  TO  FIELD  MEASUREMENTS .  78 

1.  Mercator .  78 

D.  SYNOPSIS  .  80 

V.  CONCLUSIONS  AND  RECOMMENDATIONS .  83 

APPENDIX  A.  DERIVATION  OF  ANALYTIC  SOLUTION  FOR 

THE  MUNK  PROBLEM  .  85 

LIST  OF  REFERENCES  .  91 

INITIAL  DISTRIBUTION  LIST  .  95 


viii 


LIST  OF  FIGURES 


1.  Mesh  Grid  Points  Provided  by  AWI  and  Indian  Ocean  Bathymetry  .  .  14 

2.  Initial  Conditions  for  the  Simulation  of  the  December  26,2004,  Indian 

Ocean  tsunami  [After:  [17]]  .  16 

3.  Locations  of  the  Tide  Gauges  Where  the  December  26,  2004,  Tsunami 

Waves  Were  Recorded  [After:  [22,  23,  25]]  .  18 

4.  Time  Series  of  the  Sea  Surface  Elevation  in  Paradip  -  India  [After:  [23]]  20 

5.  Results  From  De-tidal  Procedure  for  Colombo  -  Sri  Lanka  Tide  Gauge 

Records.  [After:  [25], [29], [30]] .  21 

6.  Results  From  De-tidal  Procedure  for  Hanimaadhoo  -  Maldives  Tide 

Gauge  Records.  [After:  [25], [29], [30]] .  22 

7.  Results  From  De-tidal  Procedure  for  Male  -  Maldives  Tide  Gauge  Records. 

[After:  [25], [29], [30]]  22 

8.  Results  From  De-tidal  Procedure  for  Gan  -  Maldives  Tide  Gauge  Records. 

[After:  [25], [29], [30]]  23 

9.  Results  From  De-tidal  Procedure  for  Diego  Garcia  -  UK  Tide  Gauge 

Records.  [After:  [25], [29], [30]] .  23 

10.  Results  From  De-tidal  Procedure  for  Port  Louis  -  Mauritius  Tide  Gauge 

Records.  [After:  [25], [29], [30]] .  24 

11.  Results  From  De-tidal  Procedure  for  Salalah  -  Oman  Tide  Gauge  Records. 

[After:  [25], [29], [30]]  24 

12.  Results  From  De-tidal  Procedure  for  Pointe  La  Rue  -  Seycelles  Tide 

Gauge  Records.  [After:  [25], [29], [30]] .  25 

13.  Satellites  Crosses  During  the  2004  Indian  Ocean  Tsunami  .  26 

14.  The  Mercator’s  Anchored  Location .  27 

15.  Male  Tide  Gauge  Position .  30 

16.  Comparison  of  arrival  times  -  Male  Tide  Gauge  vs.  Model .  31 


IX 


17.  Gan  Tide  Gauge  Position  .  33 

18.  Comparison  of  arrival  times  -  Gan  Tide  Gauge  vs.  Model .  34 

19.  Diego  Garcia  Tide  Gauge  Position .  35 

20.  Comparison  of  arrival  times  -  Diego  Garcia  Tide  Gauge  vs.  Model  ...  36 

21.  Port  Louis  Tide  Gauge  Position .  37 

22.  Comparison  of  arrival  times  -  Port  Louis  Tide  Gauge  vs.  Model  ....  38 

23.  Salalah  Tide  Gauge  Position .  40 

24.  Comparison  of  arrival  times  -  Salalah  Tide  Gauge  vs.  Model .  41 

25.  Pointe  La  Rue  Tide  Gauge  Position .  42 

26.  Comparison  of  arrival  times  -  Pointe  La  Rue  Tide  Gauge  vs.  Model  .  .  43 

27.  Mormugao  Tide  Gauge  Position .  44 

28.  Comparison  of  arrival  times  -  Mormugao  Tide  Gauge  vs.  Model  ....  45 

29.  Okha  Tide  Gauge  Position .  47 

30.  Comparison  of  arrival  times  -  Okha  Tide  Gauge  vs.  Model .  48 

31.  Time  Series  of  the  Sea  Surface  Elevation  in  Okha  -  India  [After:  [23]]  .  49 

32.  Paradip  Tide  Gauge  Position  .  52 

33.  Comparison  of  arrival  times  -  Paradip  Tide  Gauge  vs.  Model .  53 

34.  Model’s  Sea  Surface  Elevation  Two  Hours  After  the  Earthquake  ....  54 

35.  Model’s  Sea  Surface  Elevation  Three  Hours  After  the  Earthquake  ...  55 

36.  Vishakhapatham  Tide  Gauge  Position  .  57 

37.  Comparison  of  arrival  times  -  Vishakhapatham  Tide  Gauge  vs.  Model  .  58 

38.  Chennai  Tide  Gauge  Position .  59 

39.  Comparison  of  arrival  times  -  Chennai  Tide  Gauge  vs.  Model  .  60 

40.  Tuticorin  Tide  Gauge  Position .  62 

41.  Comparison  of  arrival  times  -  Tuticorin  Tide  Gauge  vs.  Model .  63 

42.  Kochi  Tide  Gauge  Position .  64 

43.  Comparison  of  arrival  times  -  Kochi  Tide  Gauge  vs.  Model .  65 

44.  Colombo  Tide  Gauge  Position .  67 


x 


45.  Comparison  of  arrival  times  -  Colombo  Tide  Gauge  vs.  Model .  68 

46.  Hanimaadhoo  Tide  Gauge  Position .  69 

47.  Comparison  of  Arrival  Times  Hanimaadhoo  Tide  Gauge  vs.  Model  .  .  70 

48.  Model’s  Sea  Surface  Elevation  Two  Hours  After  the  Earthquake  ....  72 

49.  Comparison  of  Jason  1  vs.  Model .  73 

50.  Model’s  Sea  Surface  Elevation  Two  Hours  and  Five  Minutes  After  the 

Earthquake .  74 

51.  Comparison  of  Topex/Poseidon  vs.  Model .  75 

52.  Model’s  Sea  Surface  Elevation  Three  Hours  and  Fifteen  Minutes  After 

the  Earthquake .  76 

53.  Comparison  of  Envisat  vs.  Model .  77 

54.  Comparison  of  arrival  times  -  Mercator’s  hshhnder  vs.  Model .  79 

55.  Energy  Propagation  of  the  Simulated  Tsunami  Waves .  81 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


LIST  OF  TABLES 


I.  Coefficients  for  the  Strong  Stability  Preserving  -  Third  Order  Runge- 

Kutta  Method .  11 

II.  Properties  of  the  Used  Mercator  Projection  .  15 

III.  Points  That  Form  the  Two  Fault  Segments  [After:  [17]] .  16 

IV.  Tsunami  Characteristics  Estimated  From  the  Indian  Tide  Gauge  Records 

[After:  [22]] .  19 

V.  Profiles  of  Sea  Surface  Height  Obtained  by  Satellites  During  the  2004 

Indian  Ocean  Tsunami  [After:  [27]] .  26 

VI.  Summary  of  Results  for  Tide  Gauge  Stations  Located  at  the  Tropical 

Indian  Ocean  and  west  of  India  regions .  50 

VII.  Summary  of  Results  for  the  Tide  Gauge  Stations  Located  at  the  South- 

Sontheast  of  India  regions .  71 


xiii 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xiv 


ACKNOWLEDGMENTS 


I  gratefully  acknowledge  the  Alfred- Wegener  Institute  for  Polar  and  Marine 
Research  (tsunami  modeling  group)  for  providing  the  real  bathymetry  data  of  the 
Indian  Ocean  region. 

I  acknowledge  the  NOAA  for  providing  satellite  altimetry  data  used  to  compare 
model  results.  I  would  also  like  to  acknowledge  the  University  of  Hawaii  Sea  Level 
Center,  Honolulu,  Hawaii,  the  Survey  of  India,  Delhi,  and  the  National  Institute  of 
Oceanography,  Goa,  India  for  providing  the  sea  level  records  of  the  2004  Indian  Ocean 
tsunami. 

I  would  like  to  thank  Alike  Cook  from  the  Oceanography  Department  of  NPS 
whose  knowledge  in  MATLAB  was  inestimable.  I  acknowledge  Rich  Pawlowicz’s 
MATLAB  codes  (m  map  and  t  tide)  which  I  used  in  order  to  convert  the  longitude 
/  latitude  coordinates  to  cartesian  and  remove  the  tidal  effect  from  the  real  sea  level 
records. 

I  would  like  to  thank  my  two  advisors,  Prof.  Frank  Giraldo  and  Prof.  Timour 
Radko,  whom  without  their  guidance  this  thesis  would  not  have  been  possible.  I  would 
also  like  to  thank  the  Naval  Postgraduate  School,  the  Oceanography  and  Mathematics 
departments  at  NPS,  and  the  Hellenic  Navy. 


xv 


THIS  PAGE  INTENTIONALLY  LEFT  BLANK 


xvi 


I. 


INTRODUCTION 


The  9.3  magnitude  Sumatra- Andaman  undersea  earthquake  occurred  at  00:58:53 
UTC  (07:58:53  local  Indonesian  time)  on  December  26,  2004.  This  earthquake  gener¬ 
ated  a  megatsunami  that  was  among  the  deadliest  disasters  in  modern  history,  killing 
more  than  226,000  people  [1]  and  is  known  as  the  2004  Indian  Ocean  tsunami. 

The  modern  science  triggered  by  this  and  other  similar  disasters  aimed  to 
predict  when  and  where  an  earthquake  would  occur  and  where  a  tsunami  might  be 
initiated.  While  the  accurate  prediction  of  an  earthquake  itself  remains  a  major 
challenge,  simulating  the  propagation  and  inundation  of  a  tsunami  wave  is  now  a 
distinct  possibility.  Knowing  where  the  faults  are  through  the  use  of  real  bathymetry 
data,  numerical  simulations  can  be  used  to  identify  the  regions  that  are  likely  to  be 
affected  by  tsunami  waves  and  establish  seismic  criteria  for  issuing  tsunami  warnings 
in  the  case  of  actual  tsunamis.  Furthermore,  by  using  the  existing  technology  of 
monitoring  seismic  activity,  the  simulation  of  the  generated  tsunami  can  give  sufficient 
warning  to  the  regions  that  are  in  the  path  of  the  wave. 

Tsunami  simulation,  as  well  as  any  other  numerical  modeling,  requires  the 
discretization  of  the  physical  space.  The  two  major  strategies  used  are  regular  and 
unstructured  discretizations  [2],  Models  that  are  based  on  regular  meshes  and  have 
been  successfully  applied  to  the  2004  Indian  Ocean  tsunami  are  as  follows:  the  MOST 
(Method  of  Splitting  Tsunami)  propagation  model  by  Titov  and  Synolakis  [3,  4], 
which  is  based  on  a  finite-difference  numerical  approximation  to  the  nonlinear  shal¬ 
low  water  wave  equations,  and  the  fully  nonlinear  and  dispersive  Boussinesq  model 
(FUNWAVE)  by  Watts  [5].  An  example  of  a  model  based  on  unstructured  grids  is 
the  TsunAWI  model  by  the  tsunami  modeling  group  of  the  Alfred- Wegener  Institute 
[2]- 

The  discontinuous  Galerkin  (DG)  method,  because  of  the  high-order  accuracy, 
geometric  flexibility  to  use  unstructured  grids,  local  conservation,  and  monotonicity 
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properties,  has  been  accepted  in  the  last  decade  by  the  geosciences  community  as  a 
critical  component  of  the  computational  geophysical  fluid  dynamics  and  specifically, 
in  the  construction  of  future  ocean  and  shallow  water  models  [6].  Schwanenberg 
and  Kongeter  [7]  first  used  the  DG  method  for  the  planar  shallow  water  equations 
followed  by  the  works  of  Li  and  Liu  [8]  as  well  as  Ainzinger  and  Dawson  [9] .  Later, 
Dupont  and  Lin  [10],  Eskilsson  and  Sherwin  [11],  Remade  et  al.  [12],  and  Kubatko 
et  al.  [13]  constructed  shallow  water  models  on  triangles  using  a  collapsed  local 
coordinate  (modal)  discontinuous  Galerkin  method.  Giraldo  et  al.  [14],  first  used  the 
DG  method  for  the  shallow  water  equations  on  the  sphere  and  later  by  construction 
on  triangular  domains  [15].  Giraldo  and  Warburton  et  al.  [6]  also  applied  the  method 
to  the  two-dimensional  oceanic  shallow  water  equations.  Motivation  for  the  choice 
to  use  unstructured  grids  is  related  to  its  ability  to  better  represent  the  continental 
coastlines  and  simulate  wave  processes  on  different  scales.  A  tsunami’s  motion  is 
characterized  by  relatively  fast  speed  and  low  amplitude  in  the  deep  ocean,  and  when 
it  reaches  shallow  waters  the  speed  decreases  and  the  amplitude  increases.  Near  the 
shoreline,  the  tsunami  steepens  and  may  break. 

Giraldo  and  Warburton  [6]  DG  formulation  used  high-order  Langrange  poly¬ 
nomials  on  the  triangle  using  nodal  sets  up  to  the  15th  order.  By  using  six  test  cases 
(three  of  which  had  analytic  solutions)  they  showed  that  the  high-order  triangular 
DG  method  exhibits  exponential  convergence. 

According  to  [16]  the  evolution  of  earthquake-generated  tsunamis  has  three 
distinctive  stages:  generation,  propagation,  and  run-up.  This  study,  as  being  part 
of  a  research  in  regards  to  the  development  of  a  triangular  discontinuous  Galerkin 
oceanic  shallow  water  model,  focuses  on  formatting  real  bathymetry  data  of  the  In¬ 
dian  Ocean  in  order  to  simulate  the  propagation  stage  of  the  Indian  Ocean  tsunami  of 
December  26,  2004,  by  using  this  DG  model.  The  initial  conditions  (i.e.,  the  tsunami’s 
generation)  are  taken  from  the  reconstruction  of  the  rupture  as  described  in  [17]  and 
based  on  [18],  while  the  run-up  stage  is  out  of  the  scope  for  this  thesis.  Real  measure- 
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merits  are  used  for  the  validation  of  this  simulation.  The  model  results  are  compared 
to  tide  gauge  data  from  several  stations  around  the  Indian  Ocean,  satellite  altimetry 
(Jason- 1,  Topex-Poseidon  and  Envisat  data),  and  field  measurements.  These  results 
show  that  the  model  gives  accurate  estimates  of  arrival  times  in  distant  locations. 

This  study  is  organized  as  follows.  Chapter  II  describes  the  spatial  discretiza¬ 
tion  of  the  governing  shallow  water  equations  and  the  time  integrator  used  in  this 
model.  Chapter  III  describes  the  real  bathymetry  data,  the  initial  conditions  used 
for  the  simulation,  and  the  real  measurement  data  used  for  comparisons.  Chapter 
IV  presents  the  results  from  the  comparison  and  Chapter  V  discusses  the  conclusions 
and  recommendations.  Appendix  A  describes  the  derivation  of  the  analytic  solution 
for  the  Munk  problem.  This  solution  is  given  in  order  for  the  validation  of  the  lateral 
diffusion  operators  in  future  studies. 
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II. 


BACKGROUND 


A.  GOVERNING  EQUATIONS 
1.  Shallow  Water  Equations 

The  oceanic  shallow  water  equations  are  a  system  of  nonlinear  partial  differen¬ 
tial  equations  which  govern  the  motion  of  a  viscous  incompressible  fluid  in  a  shallow 
depth.  The  underlying  assumption  is  that  the  depth  of  the  fluid  is  small  compared  to 
the  wave  length  of  the  disturbance.  The  Indian  Ocean  is  not  shallow  since  the  depth 
is  several  kilometers,  but  the  devastating  tsunami  on  December  26,  2004,  involved 
waves  that  were  tens  or  hundreds  of  kilometers  long,  so  the  shallow  water  approxi¬ 
mation  provides  a  reasonable  model  in  this  situation.  A  tsunami  is  the  ideal  shallow 
water  flow  problem  because  of  the  barotropic  motion  of  the  water  column. 

The  shallow  water  equations  in  conservation  form  are: 

If  +  V  ■  Flq)  =  S(q)  (II.l) 

where  q  =  (0,  (f)UT)T  are  the  conservation  variables, 

F  (?)  =  (  )  (112) 

y  4>u  <g>  u  +  \(j)2T2  —  zW  ((fm)  J 

is  the  flux  tensor,  and 

S(q)  =  -\  °  ,  T  (H.3) 

y  /  (k  x  <jm)  +  cj)'V(j)b  —  ^  +  7 <jm  J 


is  the  source  function  where: 

-  the  nabla  operator  is  defined  as  V  =  (dx,  dy)T 

-  ®  denotes  the  tensor  product  operator 

-  <j)  —  gh  is  the  geopotential  height 

(where  g  is  the  gravitational  constant  and  h  is  the  free  surface  height  of  the  fluid), 

-  (f)b  is  the  bathymetry 
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-  u  —  (u,  v)T  is  the  velocity  vector 

■  f  —  fo  +  /3(y  —  Vm)  is  the  Coriolis  parameter 

-  k  —  (0,  0, 1)T  is  the  unit  normal  vector  of  the  x  -  y  plane 

-  the  vector  r  is  the  wind  stress 

-  the  constant  7  is  the  bottom  friction 

-  the  term  X2  is  the  rank- 2  identity  matrix 

B.  TRIANGULAR  DISCONTINUOUS  GALERKIN  METHOD 

This  section  describes  the  discretization  of  the  shallow  water  equations  by  the 
discontinuous  Galerkin  method  on  triangles,  as  defined  by  Giraldo  and  Warbnrton 
[6], 

1.  Basis  Functions 

This  method  demands  the  decomposition  of  the  domain  G  into  Ne  non-overlapping 
triangular  elements  Ge  such  that 

Ne 

n=  U^e 

e=l 

and  the  introduction  of  a  nonsingular  mapping  x  =  *&(£)  which  defines  a  transforma¬ 
tion  from  the  physical  Cartesian  coordinate  system  x  =  ( x,y)T  to  the  local  reference 
coordinate  system  £  =  (£,  r])T,  defined  on  the  reference  triangle 

=  {(£,v),  -1  <  Lv  <  1,  £  +  n  <  0,}. 

The  local  element-wise  solution  q  can  be  represented  by  an  TVth  order  poly¬ 
nomial  in  ^  as 

mn 

Vn(()  =  I ><«)«at(0  (n.4) 

i= 1 

where  represents  MN  =  |(iV  +  1)(N  +  2)  interpolation  points  and  are  the 

associated  multivariate  Lagrange  polynomials. 
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The  Lagrange  polynomials  (nodal  basis  functions),  ^,(£,7 l)i  011  the  reference 
triangle  are  defined  as 


(2  i  +  1)  i}  +  j  +  1)  pQ,o 


(W) 


1  ~  V\  1  p2i+l,0 

2  J  j 


(v) 


(II  .5) 


where  -P“,/3(£)  represents  the  nth  order  Jacobi  polynomial  in  the  interval  —1  <  £  <  1, 
k  —  i  +  j(N  +  1)  +  1,  and  the  indices  vary  as  0  <  i,  j;  i  +  j  <  N,  and  k  —  1, . . . ,  M n- 
An  explicit  formula  for  the  Lagrange  basis  is 


mn 

A  (f ,  v)  =  A*Vk  (€,  v)  (11.6) 

k= 1 

where  the  indices  are  defined  as  i,  j,  k  —  1, . . . ,  M n-  By  using  the  cardinal  property 
of  the  Lagrange  polynomials 

mn 

$ij  y  ]  Aikp>k  (‘(j,  Tjj )  , 
fc=l 

where  5  is  the  Kronecker  delta  function,  it  can  be  written  as 


Aik  —  {^Pk 


(II-7) 


Recognizing  that 


bjfc  —  (fk  v j) 


(II.8) 


is  the  generalized  Vandermonde  matrix  and  using  Equations  (II. 6),  (II. 7),  and  (II. 8), 
the  Lagrange  polynomials  can  be  constructed  as  follows: 

rp 

)  =  EV‘L^K,»)  (11.9) 
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2.  Integration 

For  any  two  functions  (/)  and  (g),  the  2D  (area)  integration  X4  that  is  required 
by  the  weak  formulation  of  Galerkin  methods  is  given  by 


,  MC 

[/,g]  =  f{x)g{x)dx  =  X!  w?  |  Je(^)  |  /(&)g(&) 

J\le 

where  Mq  is  a  function  of  C  which  represents  the  order  of  the  cubature  approximation. 
For  Wi  and  £,;  the  high-order  cubature  rules  for  the  triangle  are  used. 

The  DG  method  requires  the  evaluation  of  boundary  integrals.  This  is  the 
mechanism  by  which  the  fluxes  across  element  edges  are  evaluated  and  allows  the  dis¬ 
continuous  elements  to  communicate.  According  to  [6],  the  ID  (boundary)  integration 
Tts  is  given  by 


2b  [/, g]  =  [  f(x)g(x)dx  =  I  JS(^i)  I  /(&)g(&) 

e  i=0 

where  Q  represents  the  order  of  the  quadrature  approximation.  The  use  of  Q 
gives  order  2N  —  1  accuracy. 


=  N 


3.  Semi-discrete  Equations 

The  use  of  the  discontinuous  Galerkin  discretization  in  Equation  (II.  1)  gives 
the  weak  form  of  the  DG  method 

—  Fn  ■  V  —  -0 i(x)dx  =  —  0 i(x )  n  ■  F*Ndx  (II. 10) 

where  FN  =  F(qN)  and  Sn  =  S(qN)  with  F  and  S  given  by  Equations  (II. 2)  and 
(II. 3),  respectively. 

In  the  boundary  integral  of  Equation  (II.  10)  n  is  the  outward  pointing  unit 
normal  vector  of  the  element  edge  Te  and  F*N  is  the  Rusanov  numerical  flux  as  given 
by  the  following  equation: 


n  =  v 


2  L 


Fn  (qk)  +  Fn  (q$)  -  |A|  (q%  «-  qLN 


n 


(11.11) 


where  A  =  max  (\UL\  +  \  Ufi\  +  with  UL,R  =  uL'R  ■  n  being  the  normal 

component  of  velocity  with  respect  to  the  edge  Te,  and  the  superscripts  L  and  R 


represent  the  left  and  right  sides  of  the  element  edge.  The  normal  vector  n  is  defined 
as  pointing  outward  from  left  to  right. 

Integrating  Equation  (II. 10)  by  parts  once  more  yields  the  strong  form 

JQ  ik{x)  +  V  •  Fn  —  dx  =  jf  A(x)  n  ■  (F N  —  F*N)  dx.  (11.12) 

4.  Matrix  Form  of  the  Semi-discrete  Equations 

Using  the  polynomial  approximation 

Mjv 

qn  =  Y,  ^ i 

i= 1 

the  semi-discrete  system  can  be  written  as 


f  4’i4’j  dec  +  /  ipi'V'ipj  dx  ■  F j  —  f  dxSj  =  f  ipiipjndx  ■  (F  —  F*)r 
*/ ilg  (/f  « lip  « Og  JTe 

(11.13) 

Next,  the  elemental  matrices  Mf,  are  dehned  as  follows: 

'  LJ 

Mij  =  Jn  jdx ,  M\j  =  jf  jTidx ,  D'j  =  ifiiViftjdx, 

which  makes  it  possible  to  write  Equation  (11.13)  in  the  matrix  form: 


M, 


v 


dt 


D%)  Fj  -  MijSj  —  (M(.)  (F  -  F 


(11.14) 


where  the  superscript  e  denotes  an  element-based  evaluation  and  s  denotes  a  side- 
based  (or  edge-based)  evaluation. 

The  elemental  matrices  can  be  written,  in  terms  of  the  computational  variables 
(£),  as  follows: 

Mtj  =  \r\  L  ^3di  =  \r\Mi3 

d  '  fag 


r 


At  =  \Je\(Dt,G  +  Dlr,l)i  +  \r\(D^l  +  Dl-nl)] 


u  >y) 


Ms-  = 

y 


Js 


ipi^jndi  =  \Js\M*Jnsxi  +  nsj 
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where  the  metric  terms  are  defined  as  \Je\  =  2Ae  and  |JS|  =  2AS.  The  Ae  and  As 
denote  the  area  of  element  e  and  the  length  of  edge  s  respectively.  Note  that  the 
gradient  operator  in  these  integrals,  V^,  is  defined  in  the  local  reference  coordinate 
system  £,  and  Vte  and  Te  denote  the  area  and  boundary  domains  in  the  computational 
space,  namely  the  bounds  of  integration  for  the  master  element.  Thus,  Equation 
(11.14)  can  be  written  in  the  following  way: 


+  \r\ (4e  +  DH)  r,  +  \r\ (d^i + ak)  -  \je\ M»si 

(11.15) 


=  K(r-r),  +  nl(g‘ -g‘), 


where  Fe  =  fei  +  gej. 

By  dividing  Equation  (11.15)  by  | ,Je\  and  left  multiplying  by  M the  following 
is  obtained 


ggj 

dt 


+  (die + Die)  r, + (Die + DU 


v.  Si 


Sf  = 


\Je 


[Mi 


v 


nUr-Dj  +  Ktf-glj 


where  the  matrices  are  defined  as 


(11.16) 


where 


D%  =  M^Dl,,  E%  =  M^Dl,  Ml  =  Mg1  M‘k. 

MC  Mq 

Mij  =  ^kAkfpjk,  wkAk^jk, 

k= 1  k= 1 


(11.17) 


(11.18) 


r,S  V?  ,  ™  ^  ,  Si’ik 

Dii  -  ID  -  y  — .  (11,19) 

The  Mq  and  Mq  denote  the  number  of  cubature  (two-dimensional)  and  quadrature 
(one-dimensional)  integration  points  required  to  achieve  order  2 N  —  1  accuracy,  and 
^ ik  represents  the  function  ^  at  the  i  =  1  ,...,Mjv  interpolation  points  evaluated  at 
the  integration  point  k. 
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Since  the  mass  matrix  is  constant  (i.e. ,  not  a  function  of  x )  then  the  use  of 
Equations  (11.18)  and  (11.19)  gives  the  following: 

Mij  =  wkAkil>jki  wk$ik 

k=l  k= 1  k= 1 

where 

A  =  M^ipk 

is  the  basis  function  premultiplied  by  the  inverse  mass  matrix.  Absorbing  the  mass 
matrix  within  the  test  function  completely  eliminates  the  mass  matrix  from  Equa¬ 
tion  (11.16)  without  making  any  approximations. 


5.  Time  Integrator 

In  order  to  advance  the  solution  in  time  while  retaining  high-order  accuracy, 
the  strong  stability  preserving  (SSP)  Runge-Kutta  third-order  RK3  method  (see  [6]) 
is  used.  The  semi-discrete  (in  space)  equations  are  written  as  follows: 


dq 

dt 


=  S{q). 


The  SSP  temporal  discretization  of  this  vector  equation  is 


fork  =  1, ...,  3  : 

qk  =  akqn +  akqk~1  +  pkAtS{qk~1) 


where  q°  =  qn ,  q  ‘  =  qn+1,  and  the  coefficients  a  and  f3  are  given  in  Table  I. 


k=l 

k=2 

k=3 

«0 

1 

3/4 

1/3 

on 

0 

1/4 

2/3 

p 

1 

1/4 

2/3 

Table  I.  Coefficients  for  the  Strong  Stability  Preserving  -  Third  Order  Runge-Kutta 
Method. 
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III.  SIMULATION  OF  THE  INDIAN  OCEAN 

TSUNAMI 

A.  MODEL  DESCRIPTION 

1.  Boundary  Conditions 

As  was  mentioned  above,  this  thesis  study  focuses  on  the  simulation  of  the 
propagation  stage  of  the  Indian  Ocean  tsunami  while  the  inundation  modeling  is  out 
of  its  scope.  For  this  reason,  the  shoreline  is  being  treated  as  a  fixed-wall  boundary. 
These  no-flux  boundary  conditions  are  enforced  at  all  land  boundaries  by  virtue  of 
the  statement 

n  ■  u  =  0  (III.  1) 

which  eliminates  the  normal  component  of  the  velocity  to  the  no-flux  boundary  with¬ 
out  altering  the  tangential  component  (also  known  as  free  slip  boundary  conditions). 
Furthermore,  the  points  near  the  coastlines  with  depth  less  than  11  meters  are  con¬ 
sidered  to  have  11  meters  depth. 

Since  the  model  does  not  currently  have  inundation  algorithms,  it  only  makes 
sense  to  compare  observations  with  the  model  results  using  arrival  times  and  heights 
of  the  first  waves.  This  way  the  model  results  will  not  be  affected  by  the  errors 
associated  with  representation  of  boundary  conditions  along  the  coastlines. 

2.  Tides  /  Coriolis  /  Viscosity 

This  study  is  based  on  two  experiments.  Each  experiment  had  ten  hours  of 
simulations  with  time  steps  of  1.5  seconds.  During  the  simulations,  tidal  forcing 
and  the  Coriolis  effect  were  not  included.  In  the  first  experiment,  the  value  of  the 
horizontal  viscosity  coefficient  was  v  =  1000,  while  in  the  second  the  value  was  zero. 
Both  simulations  yielded  the  same  results.  This  can  be  explained  by  using  scale 
analysis  for  the  shallow  water  momentum  equations,  which  indicates  that  the  relation 
of  the  horizontal  viscosity  to  the  rest  of  the  forces  is  10-8,  on  the  order  of  machine 
single  precision. 
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B.  REAL  BATHYMETRY  DATA 


For  this  study,  real  bathymetry  data  of  the  Indian  Ocean  was  provided  by  the 
Alfred- Wegener  Institute  for  Polar  and  Marine  Research  (tsunami  modeling  group). 
The  unstructured  mesh  was  produced  with  the  mesh  generator  TRIANGLE  by  Jonathan 
Shewchuk  [20].  As  described  in  [2],  the  resulting  meshes  were  smoothed  to  improve 
the  overall  mesh  quality.  The  coarsest  resolution  of  the  mesh  is  14  kilometers  in  the 
deep  ocean,  500  meters  in  the  coastal  areas  and  in  areas  with  high  bathymetry  gra¬ 
dients  (canyons),  and  less  than  100  meters  in  the  northern  tip  of  Sumatra.  The  total 
number  of  elements  used  was  130, 445  created  by  66,  715  grid  points.  Figure  1  presents 
the  grid  points  used  for  this  simulation  study  and  the  Indian  Ocean  bathymetry  in 
meters. 


Longitude 


4000 

3000 

2000 

1000 

0 

-1000 

-2000 

-3000 

-4000 

-5000 


Figure  1.  Mesh  Grid  Points  Provided  by  AWI  and  Indian  Ocean  Bathymetry 


The  conversion  of  the  above  real  bathymetry  data  from  latitude/longitude  co¬ 
ordinates  into  Cartesian  coordinates  (x,y)  was  made  by  using  the  MATLAB  function 
m_U2xy  as  described  in  [21].  The  used  projection  was  a  Mercator  projection  with 
the  properties  given  in  Table  II. 
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Min  Longitude 

Max  Longitude 

Min  Latitude 

Max  Latitude 

030°U 

110°U 

35°  S 

35°N 

Table  II.  Properties  of  the  Used  Mercator  Projection 


C.  TSUNAMI  SOURCE  MODEL 

The  main  cause  of  tsunami  waves  is  the  abruptly  motion  of  converging  or 
destructive  plate  boundaries  which  vertically  displace  the  overlying  water.  The  Indian 
Ocean  tsunami  was  generated  by  the  static  sea  floor  uplift  caused  by  an  abrupt  slip 
of  the  India/Burma  plate  interface.  This  permanent,  vertical  sea  floor  displacement 
can  be  computed  using  the  static  dislocation  formulae  from  [18].  Inputs  for  this 
computation  are  the  fault  plane  location,  its  depth,  strike,  dip,  slip,  length,  width, 
seismic  moment,  and  rigidity.  The  initial  conditions  for  the  tsunami  simulation  were 
taken  from  the  reconstruction  of  the  rupture  as  described  in  [17].  According  to  this 
study,  the  fault  extent  constrained  by  observed  tsunami  arrival  time  to  the  northwest, 
east  and  south  of  the  slip  zone  indicates  a  fault  zone  of  approximately  1000  kilometers 
by  200  kilometers.  The  epicenter  location  lies  on  the  southern  end  of  the  fault  zone. 
To  accommodate  trench  curvature,  this  fault  plane  was  broken  into  two  segments  as 
depicted  in  Figure  2.  The  points  that  form  these  planes  and  used  for  this  study  are 
presented  in  Table  III. 
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Figure  2.  Initial  Conditions  for  the  Simulation  of  the  December  26,2004,  Indian  Ocean 
tsunami  [After:  [17]] 


Fault  Segment 

Corner 

Latitude 

Longitude 

Sea  surface 
elevation 

Northern  Fault  Segment 

NW 

11. 78°  AT 

092.12°E 

+5.07 

nr 

Northern  Fault  Segment 

sw 

5.6°N 

093.22°E 

+5.07 

m 

Northern  Fault  Segment 

NE 

12.05  °N 

094.02°E 

-4.75 

nr 

Northern  Fault  Segment 

SE 

6.0°N 

094.97°E 

-4.75 

m 

Southern  Fault  Segment 

NW 

5.33°N 

093.25°E 

+5.07 

nr 

Southern  Fault  Segment 

SW 

2.97°  N 

094.35°E 

+5.07 

m 

Southern  Fault  Segment 

NE 

6.0°N 

094.97°E 

-4.75 

m 

Southern  Fault  Segment 

SE 

3.88°N 

095.97°E 

-4.75 

m 

Table  III.  Points  That  Form  the  Two  Fault  Segments  [After:  [17]] 
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D.  REAL  MEASUREMENT  DATA 
1.  Tide  Gauge  Records 

Tide  gauge  records  are  often  used  to  provide  valuable  information  on  the 
tsunami  arrival  times  and  the  changes  in  sea  level  due  to  tsunamis.  After  re¬ 
moving  the  tidal  effect  these  data  provide  the  changes  in  water  level  due  to  tsunamis 
alone. 

The  Indian  Ocean  tsunami  was  recorded  on  all  the  tide  gauges  located  in  the 
Indian  Ocean.  A  comprehensive  overview  of  these  stations  and  an  analysis  of  the 
tsunami  records  is  published  in  [22],  For  this  thesis  study  tide  gauge  records  from 
fifteen  stations  in  the  Indian  Ocean  were  used.  The  locations  of  these  stations  are 
shown  in  Figure  3.  In  order  to  increase  the  accuracy  of  the  results,  each  station’s 
location  was  interpolated  using  the  three  grid  points  of  the  triangular  element  in 
which  the  station  belongs.  The  tsunami  characteristics  of  each  station’s  records  are 
shown  in  Table  IV. 
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Figure  3.  Locations  of  the  Tide  Gauges  Where  the  December  26,  2004,  Tsunami 
Waves  Were  Recorded  [After:  [22,  23,  25]] 
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No 

Station 

Coordi¬ 

nates 

Sampling 

interval 

(min) 

First 

Arrival 
time  (UTC) 

wave 

Travel 

time 

1 

Paradip 

20.26 °N 
086.70 °E 

6 

03  :  27 

2 

hrs  28  min 

2 

Vishakhapatham 

17.65  °N 
083.28 °E 

5 

03  :  35 

2 

hrs  36  min 

3 

Chennai 

13.10° N 
080.32 °E 

5 

03  :  33 

2 

hrs  34  min 

4 

Tnticorin 

08.75 °N 
078.20 °E 

6 

04  :  24 

3 

hrs  25  min 

5 

Kochi  (Cochin) 

09.97 °N 
076.27 °E 

6 

05  :  41 

4 

hrs  42  min 

6 

Mormugao 

15.42°A 

073.80°A 

5 

06  :  53 

5 

hrs  54  min 

7 

Okha 

22.50 °N 
069.10°A 

6 

09  :  03 

8 

hrs  04  min 

8 

Colombo 

06.93°A 

079.83°A 

2 

03  :  49 

2 

hrs  50  min 

9 

Hanimaadhoo 

06. 77°  A 
073.18°F 

2 

04  :  30 

3 

hrs  31  min 

10 

Male 

04.18°A 

073.52°A 

4 

04  :  14 

3 

hrs  15  min 

11 

Gan 

00.68°A 
073. 17°  E 

4 

04  :  16 

3 

hrs  17  min 

12 

Diego  Garcia 

07.30 °S 
072.38 °E 

6 

04  :  45 

3 

hrs  46  min 

13 

Port  Lonis 

20.15  °A 
057.50°A 

2 

07  :  46 

6 

hrs  47  min 

14 

Salalah 

17.00°V 

054.00°A 

4 

08  :  08 

7 

hrs  09  min 

15 

Pointe  La  Rue 

04.68°A 

055.53°A 

4 

08  :  16 

7 

hrs  17  min 

Table  IV.  Tsunami  Characteristics  Estimated  From  the  Indian  Tide  Gauge  Records 
[After:  [22]] 
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The  records  of  the  first  seven  stations  of  Table  IV  (these  stations  are  depicted 
with  their  names  in  black  color  in  Figure  3)  were  obtained  from  the  National  Institute 
of  Oceanography  (NIO),  Goa,  India  [23].  These  tide  gauges  are  part  of  the  network 
maintained  by  the  Survey  of  India  (SOI)  agency  along  the  coast  of  India.  According 
to  [24],  SOI  tide  gauges  were  either  mechanical  float-type  analog  gauges  or  pressure¬ 
sensor  gauges.  The  analog  records  from  Vishakhapatham,  Chennai,  and  Mormugao 
were  digitized  by  SOI  at  an  interval  of  5  minutes,  while  those  from  Okha  at  an  interval 
of  6  minutes.  The  pressure-sensor  gauges  from  Paradip,  Tuticorin,  and  Kochi  had 
sampling  intervals  of  6  minutes.  All  these  are  high-quality,  de-tided  records.  Some 
extensive  gaps  in  the  data  were  due  to  instruments  and  transmission  problems.  A 
typical  example  of  these  records  is  presented  in  Figure  4. 


Paradip  -  Observation  Data  -  Pos.:  20.26N  ;  086. 70E 


Figure  4.  Time  Series  of  the  Sea  Surface  Elevation  in  Paradip  -  India  [After:  [23]] 
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The  records  of  the  next  eight  stations  of  Table  IV  (these  stations  are  depicted 
with  their  names  in  red  color  in  Figure  3)  were  obtained  from  the  University  of 
Hawaii’s  Sea  Level  Center  database  (Honolulu)  [25].  These  are  digital  Global  Sea 
Level  Observing  System  (GLOSS)  stations  with  2  minutes  sampling  intervals  for 
Colombo,  Hanimaadhoo,  and  Port  Louis;  4  minutes  for  Male,  Gan,  Salalah,  and 
Pointe  La  Rue;  and  6  minutes  for  Diego  Garcia.  Most  of  these  stations  are  located  on 
isolated  islands  and  therefore  the  records  were  not  significantly  affected  by  reflections 
[22],  These  records  were  real-time  sea  level  measurements  (that  included  tides)  and 
were  de-tided  by  using  the  t-tide  program  as  described  in  [29]  and  [30].  Two  problems 
during  this  procedure  were  the  extensive  gaps  in  the  data  due  to  instruments  and 
transmission  problems  and  the  limited  total  record  time  (only  one  week).  The  results 
are  presented  in  Figures  5  through  12.  In  all  the  above  cases,  except  for  Diego  Garcia, 
the  tidal  effect  was  significantly  reduced. 


Colombo  -  Observation  Data  -  Pos.:  06.93N  ;  079.83E  Colombo  -  De-Tidal  Data  -  Pos.:  06.93N  ;  079.83E 


Figure  5.  Results  From  De-tidal  Procedure  for  Colombo  -  Sri  Lanka  Tide  Gauge 
Records.  [After:  [25], [29], [30]] 
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Hanimaadhoo  -  Observation  Data  -  Pos.:  06.77N  ;  073.1 8E 


Hanimaadhoo  -  De-Tidal  Data  -  Pos.:  06.77N  ;  073.1 8E 


Figure  6.  Results  From  De-tidal  Procedure  for  Hanimaadhoo  -  Maldives  Tide  Gauge 
Records.  [After:  [25], [29], [30]] 


Male  -  Observation  Data  -  Pos.:  04.1 8N  ;  073.52E 


Male  (De-Tidal)  -  Pos.:  04.1 8N  ;  073.52E 


Figure  7.  Results  From  De-tidal  Procedure  for  Male  -  Maldives  Tide  Gauge  Records. 
[After:  [25], [29], [30]] 
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Gan  -  Observation  Data  -  Pos.:  00.68S  ;  073.1 7E 


Gan  -  De-Tidal  -  Pos.:  00.68S  ;  073.1 7E 
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Figure  8.  Results  From  De-tidal  Procedure  for  Gan  -  Maldives  Tide  Gauge  Records. 
[After:  [25], [29], [30]] 


Diego  Garcia  -  Observation  Data  -  Pos.:  07.30S  ;  072.38E  Diego  Garcia  -  De-Tidal  Data  -  Pos.:  07.30S  ;  072.38E 


Figure  9.  Results  From  De-tidal  Procedure  for  Diego  Garcia  -  UK  Tide  Gauge 
Records.  [After:  [25], [29], [30]] 
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Port  Louis  -  Observation  Data  -  Pos.:  20.1 5S  ;  057.50E 


Port  Louis  -  De-Tidal  Data  -  Pos.:  20.1 5S  ;  057.50E 
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Figure  10.  Results  From  De-tidal  Procedure  for  Port  Louis  -  Mauritius  Tide  Gauge 
Records.  [After:  [25], [29], [30]] 


Salalah  -  Observation  Data  -  Pos.:  17.00N  ;  054.00E  Salalah  -  De-Tidal  Data  -  Pos.:  17.00N  ;  054.00E 


Figure  11.  Results  From  De-tidal  Procedure  for  Salalah  -  Oman  Tide  Gauge  Records. 
[After:  [25], [29], [30]] 
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Pointe  La  Rue  -  Observation  Data  -  Pos.:  04.68S  ;  055.53E  Pointe  La  Rue  -  De-Tidal  Data  -  Pos.:  04.68S  ;  055.53E 


Figure  12.  Results  From  De-tidal  Procedure  for  Pointe  La  Rue  -  Seycellcs  Tide  Gauge 
Records.  [After:  [25], [29], [30]] 


2.  Satellite  Altimetry 

Radar  altimeters  on  board  the  Jason-1,  TOPEX/Poseidon,  and  Envisat  satel¬ 
lites  obtained  profiles  of  sea  surface  height  on  trajectories  across  the  Indian  Ocean 
between  two  and  fonr  hours  after  the  Sumatra  earthquake  on  December  26,  2004. 

Jason-1  crossed  the  equator  at  2:55  UTC,  approximately  two  hours  after  the 
earthquake  [26].  The  data  were  received  hours  to  days  after  ’’real  time”  so  it  was 
too  late  to  detect  and  warn  about  the  upcoming  tsunami,  but  it  can  be  used  to 
validate  new  models.  The  data  used  for  this  study  were  obtained  from  NOAA  [27] 
and  originated  from  the  satellite  crossings  of  Figure  13  and  Table  V. 
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Figure  13.  Satellites  Crosses  During  the  2004  Indian  Ocean  Tsunami 


No 

Satellite 

Time  after  the  Earthquake 

Color  in  Figure  13 

1 

Jasonl 

02  hrs  00  min 

Red 

2 

TOPEX  /  Poseidon 

02  hrs  05  min 

Black 

3 

Envisat 

03  hrs  15  min 

Green 

Table  V.  Profiles  of  Sea  Surface  Height  Obtained  by  Satellites  During  the  2004  Indian 
Ocean  Tsunami  [After:  [27]] 
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3.  Field  Measurements 

During  the  Indian  Ocean  tsunami  a  Belgian  yacht,  the  Mercator,  was  anchored 
about  1.6  kilometers  off  the  Phuket  coast  (Thailand)  at  07.715°IV,  098.28°i7  (Figure 
14).  The  yacht’s  depth  gauge  was  operational  and  measured  changing  wave  heights 
during  the  tsunami  [28] .  Based  on  these  measurements,  the  three  main  tsunami  waves 
had  trough-to-crest  wave  heights  of  6.6,  2.2  and  5.5  meters.  The  first  wave  (trough) 
struck  the  yacht’s  location  at  02:38  UTC  (1  hour  and  39  minutes). 


Mercator  -  Pos.:  07.71 5S  098. 28E 


Figure  14.  The  Mercator’s  Anchored  Location 
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IV. 


SIMULATIONS  RESULTS 


A.  COMPARISON  TO  TIDE  GAUGE  RECORDS 

The  simulation  results  are  separated  into  the  following  two  groups.  The  first 
group  represents  simulations  of  sea  surface  height  measurements  at  the  tropical  Indian 
Ocean  (Male,  Gan,  Diego  Garcia,  Port  Louis,  and  Pointe  La  Rue),  Oman  (Salalah) 
and  west  India  tide  gauge  stations  (Okha  and  Mormugao).  The  second  group  rep¬ 
resents  simulations  of  sea  surface  height  measurements  at  the  east  Indian  tide  gauge 
stations  (Paradip,  Vishakhapatham,  Chennai,  Tuticorin,  Kochi,  and  Colombo)  as 
well  as  at  Hanimaadhoo  (Maldives). 

1.  Tide  Gauge  Stations  Located  at  the  Tropical  Indian 
Ocean  and  West  of  India  Regions 

All  these  stations,  because  of  their  position,  were  directly  affected  by  the 
tsunami  waves  emanating  from  the  source  area. 

a.  Male  Station  (Maldives) 

This  station  was  located  at  04.18°iV,  073.52°i?  and  was  approximately 
1, 180  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  15).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  4  minutes  intervals. 

Figure  16  illustrates  the  agreement  between  the  tide  gauge  record  and 
the  model  results  about  the  arrival  time  (3  hours  and  15  minutes  after  the  earthquake) 
and  the  sign  (positive  -  wave  crest)  of  the  first  wave.  There  is  a  difference  (100  instead 
of  125  centimeters)  on  its  maximum  amplitude. 
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Latitude 


Male  -  Pos.:  04.1 8N  073. 52E 


Figure  15.  Male  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 


Tide  gauge  record  Male  vs.  Model 


Figure  16.  Comparison  of  arrival  times  -  Male  Tide  Gauge  vs.  Model 
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b.  Gan  Station  (Maldives) 

This  station  was  located  at  00.68°5',  073.17°i?  and  was  approximately 


1,260  nautical  miles  from  the  west  side  of  the  south  fault  segment  (Figure  17).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  4  minutes  intervals. 

From  Figure  18  it  is  obvious  that  the  model  correctly  predicts  the  ar¬ 
rival  time  (3  hours  and  17  minutes  after  the  earthquake),  the  sign  (positive  -  wave 
crest),  and  almost  the  amplitude  of  the  first  wave  (60  instead  of  75  centimeters). 

c.  Diego  Garcia  Station  (UK) 

This  station  was  located  at  07.30°S',  072.38°f7  and  was  approximately 
1,460  nautical  miles  from  the  west  side  of  the  south  fault  segment  (Figure  19).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  6  minutes  intervals. 

Figure  20  documents  the  agreement  between  the  tide  gauge  record  and 
the  model  results  about  the  arrival  time  (3  hours  and  46  minutes  after  the  earth¬ 
quake)  and  the  sign  (positive  -  wave  crest)  of  the  first  wave,  although  there  is  a  small 
difference  (35  instead  of  42  centimeters)  on  its  maximum  amplitude. 

d.  Port  Louis  Station  (Mauritius) 

This  station  was  located  at  20.15°S',  057.50°f7  and  was  approximately 
2,600  nautical  miles  from  the  west  side  of  the  south  fault  segment  (Figure  21).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  2  minutes  intervals,  ceased  to  operate 
for  1  hour  [22],  and  the  sign  of  the  first  wave  was  unknown. 

Figure  22  depicts  the  agreement  between  the  tide  gauge  record  and  the 
model  on  the  arrival  time  (6  hours  and  47  minutes  after  the  earthquake)  of  the  first 
wave.  According  to  the  model,  the  first  wave  was  positive  (wave  crest). 


32 


Latitude 


Gan  -  Pos.:  00.68S  073.1 7E 


Figure  17.  Gan  Tide  Gauge  Position 
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Arrival  Time 

Tide  gauge  record  Gan  vs.  Model 


Figure  18.  Comparison  of  arrival  times  -  Gan  Tide  Gauge  vs.  Model 
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Latitude 


Diego  Garcia  -  Pos.:  07. 30S  072. 38E 


Figure  19.  Diego  Garcia  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Diego  Garcia  vs.  Model 


Figure  20.  Comparison  of  arrival  times  -  Diego  Garcia  Tide  Gauge  vs.  Model 
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Latitude 


Port  Louis  -  Pos.:  20.1 5S  057. 50E 


Figure  21.  Port  Louis  Tide  Gauge  Position 
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Arrival  Time 

Tide  gauge  record  Port  Louis  vs.  Model 


Figure  22.  Comparison  of  arrival  times  -  Port  Louis  Tide  Gauge  vs.  Model 
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e.  Salalah  Station  (Oman) 

This  station  was  located  at  17.00°1V,  054.00°i7  and  was  approximately 
2,550  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  23).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  4  minutes  intervals. 

Figure  24  depicts  the  general  agreement  between  the  tide  gauge  records 
and  the  model. 


/.  Pointe  La  Rue  Station  (Seychelles) 

This  station  was  located  at  04. 68°S,  055.53°i7  and  was  approximately 
2,350  nautical  miles  from  the  west  side  of  the  south  fault  segment  (Figure  25).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  4  minutes  intervals. 

From  Figure  26  it  is  obvious  that  the  model  correctly  predicts  the  ar¬ 
rival  time  (7  hours  and  04  minutes  after  the  earthquake),  the  sign  (positive  -  wave 
crest),  and  nearly  the  amplitude  of  the  first  wave  (50  instead  of  80  centimeters). 

g.  Mormugao  Station  (India) 

This  station  was  located  at  15.42°iV  —  073.80°i7  and  was  approximately 
1,300  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  27).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  5  minutes  intervals. 

Figure  28  depicts  the  agreement  on  the  arrival  time  (5  hours  and  54  min 
after  the  earthquake)  and  the  sign  (positive  -  wave  crest)  of  the  first  wave  between  the 
tide  gauge  record  and  the  model.  Also,  the  maximum  amplitude  of  the  first  arrival 
wave  (almost  55  centimeters)  was  accurately  predicted. 
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Latitude 


Salalah  -  Pos.:  17.00N  054.00E 


Figure  23.  Salalah  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 


Tide  gauge  record  Salalah  vs.  Model 


Figure  24.  Comparison  of  arrival  times  -  Salalah  Tide  Gauge  vs.  Model 
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Latitude 


Pointe  La  Rue  -  Pos.:  04.68S  055. 53E 


Figure  25.  Pointe  La  Rue  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Pointe  La  Rue  vs.  Model 


Figure  26.  Comparison  of  arrival  times  -  Pointe  La  Rue  Tide  Gauge  vs.  Model 


43 


Latitude 


Mormugao  -  Pos.:  15.42N  073. 80E 


Figure  27.  Mormugao  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Mormugao  vs.  Model 


Figure  28.  Comparison  of  arrival  times  -  Mormugao  Tide  Gauge  vs.  Model 
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h.  Okha  Station  (India) 

This  station  was  located  at  22.50°1V,  069.10°Fi  and  was  approximately 
1,900  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  29).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  6  minutes  intervals. 

The  negative  values  of  the  observations  in  Figure  30  can  be  attributed 
to  the  tidal  effect,  which  was  not  completely  removed  (see  Figure  31).  Taking  into 
account  this  deficiency  of  data,  the  model  correctly  predicts  the  sign  (positive  -  wave 
crest),  arrival  time  (8  hours  and  04  minutes)  and  amplitude  of  the  first  wave  (8  cen¬ 
timeters)  . 
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Latitude 


Okha-  Pos.:  22.50N  069. 10E 


Figure  29.  Okha  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Okha  vs.  Model 


Figure  30.  Comparison  of  arrival  times  -  Okha  Tide  Gauge  vs.  Model 
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Sea  surface  elevation  [cm] 


Okha  -  Observation  Data  -  Pos.:  22.50N  ;  069. 10E 


Figure  31.  Time  Series  of  the  Sea  Surface  Elevation  in  Okha  -  India  [After:  [23]] 
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A  summary  of  the  simulation  results  for  the  tide  gauge  stations  located 
at  tropical  Indian  Ocean  and  west  of  India  regions  are  presented  in  Table  VI. 


No 

Tide  Gauge 
Station/ 

Field  measurement 

Coordinates 

Arrival  Time 
Difference 
(min) 

Max  wave’s  height 
Relative  error 

1 

Male 

04.18°A 

073.52°F 

0 

-0.2 

2 

Gan 

00.68°S 

073.17°F 

0 

-0.2 

3 

Diego  Garcia 

07.30°A 

072.38°F 

0 

-0.16 

4 

Port  Louis 

20.15°S' 

057.50°F 

0 

(?) 

5 

Salalah 

17.00°iV 

054.00°F 

0 

-0.1 

6 

Pointe  La  Rue 

04.68°S 

055.53°F 

0 

-0.35 

7 

Mormugao 

15.42UA 

073.80°F 

0 

-0.09 

8 

Okha 

22.50°V 

0 

0 

069.10°F 

Table  VI.  Summary  of  Results  for  Tide  Gauge  Stations  Located  at  the  Tropical  Indian 
Ocean  and  west  of  India  regions 


Comments: 

(1)  The  tide  gauge  at  Port  Louis  ceased  to  operate  for  one  hour  [22] 

(2)  The  relative  error  is  evaluated  as  follows:  Error  = 
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2.  Tide  Gauge  Stations  Located  South-southeast  of  In¬ 
dia 

The  tsunami  waves  that  arrived  at  the  region  of  these  stations  emanating  from 
the  source  (traveled  in  a  parallel  direction  to  the  earthquake  rupture)  and  to  the  in¬ 
fluence  of  tsunami  waves  reflected  from  the  southwestern  coasts  of  the  Indonesian 
Islands. 


a.  Paradip  Station  (India) 

This  station  was  located  at  20.26°1V,  086.70°77  and  was  approximately 
600  nautical  miles  from  the  north  side  of  the  north  fault  segment  (Figure  32).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  6  minutes  intervals. 

Figure  33  illustrates  the  agreement  with  respect  to  the  arrival  time 
(2  hours  and  28  minutes  after  the  earthquake)  and  the  positive  sign  (wave  crest) 
of  the  first  wave  between  the  tide  gauge  record  and  the  model.  However,  there  is 
a  significant  difference  (25  instead  of  125  centimeters)  in  the  maximum  amplitudes 
between  the  observed  height  and  those  predicted  by  the  model.  The  waves  that  arrived 
at  this  region  represent  a  combination  of  waves  emanated  from  the  source  (traveling 
parallel  to  the  earthquake  rupture)  and  waves  reflected  from  the  southwestern  coasts 
of  the  Indonesian  Islands.  Figures  34  and  35  are  snapshots  of  the  model’s  sea  surface 
elevation  two  and  three  hours  after  the  earthquake,  respectively.  These  figures  depict 
the  waves  generated  by  reflections  from  the  southwestern  coasts  of  the  Indonesian 
Islands.  Since  this  model  does  not  have  inundation  algorithms,  it  cannot  be  expected 
to  accurately  reproduce  the  maximum  height  of  the  tide  gauge  records. 
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Latitude 


Paradip  -  Pos.:  20.26N  086.70E 


Figure  32.  Paradip  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Paradip  vs.  Model 


Figure  33.  Comparison  of  arrival  times  -  Paradip  Tide  Gauge  vs.  Model 
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Figure  34.  Model’s  Sea  Surface  Elevation  Two  Hours  After  the  Earthquake 
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Figure  35.  Model’s  Sea  Surface  Elevation  Three  Hours  After  the  Earthquake 
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b.  Vishakhapatham  Station  (India) 

This  station  was  located  at  17.65°1V,  083.28°i7  and  was  approximately 


700  nautical  miles  from  the  north  side  of  the  north  fault  segment  (Figure  36).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  5  minutes  intervals. 

From  Figure  37  it  is  obvious  that  the  model  correctly  predicts  the  pos¬ 
itive  sign  (wave  crest)  of  the  first  wave  but  an  earlier  arrival  time  (2  hours  and  26 
minutes  instead  of  2  hours  and  36  minutes  after  the  earthquake),  and  underestimates 
the  amplitude  of  the  first  wave  (75  instead  of  125  centimeters).  As  in  Paradip  case, 
this  region  was  affected  by  a  combination  of  waves  from  the  source  and  waves  reflected 
from  the  southwestern  coasts  of  the  Indonesian  Islands. 

c.  Chennai  Station  (India) 

This  station  was  located  at  13.10°1V,  080.32°i7  and  was  approximately 
750  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  38).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  5  minutes  intervals. 

Figure  39  illustrates  the  agreement  on  the  arrival  time  (2  hours  and 
34  minutes  after  the  earthquake)  and  the  positive  sign  (wave  crest)  of  the  first  wave 
between  the  tide  gauge  record  and  the  model.  The  model,  however,  overestimates  the 
maximum  amplitude  (125  instead  of  70  centimeters). 
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Vishakhapatham  -  Pos.:  17.65N  083. 28E 


Figure  36.  Vishakhapatham  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 


Tide  gauge  record  Vishakhapatham  vs.  Model 


Figure  37.  Comparison  of  arrival  times  -  Vishakhapatham  Tide  Gauge  vs.  Model 


Latitude 


Chennai  -  Pos.:  13. ION  080.32E 


Figure  38.  Chennai  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Chennai  vs.  Model 


Figure  39.  Comparison  of  arrival  times  -  Chennai  Tide  Gauge  vs.  Model 
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d.  Tuticorin  Station  (India) 

This  station  was  located  at  08.75°1V,  078.20°i7  and  was  approximately 


900  nautical  miles  from  the  west  side  of  the  north  fault  segment  Figure(40).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  6  minutes  intervals. 

From  Figure  41  it  is  obvious  that  the  model  correctly  predicts  the  pos¬ 
itive  sign  (wave  crest)  of  the  first  wave  but  an  earlier  arrival  time  (3  hours  and  15 
minutes  instead  of  3  hours  and  25  minutes  after  the  earthquake),  and  underestimates 
the  amplitude  of  the  first  wave  (50  instead  of  90  centimeters).  The  earlier  arrival  time 
of  this  case  can  be  explained  by  the  fact  that  the  arrival  waves,  during  their  path 
from  the  initialization  to  the  position  of  the  gauge,  had  to  pass  through  areas  with 
depths  less  than  11  meters.  Due  to  the  boundary  conditions  used  in  the  simulations 
(all  positions  with  depth  less  than  11  meters  were  treated  as  having  a  depth  of  11 
meters  in  order  to  avoid  negative  water  depth,  this  is  due  to  the  absence  of  inundation 
algorithms  in  the  current  version  of  the  model)  the  model  predicts  earlier  arrival  times. 

e.  Kochi  Station  (India) 

This  station  was  located  at  09. 97° IV,  076.27°i7  and  was  approximately 
1,050  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  42).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  6  minutes  intervals. 

Figure  43  shows  that  the  model  correctly  predicts  the  positive  sign 
(wave  crest)  of  the  first  wave  but  an  earlier  arrival  time  (4  hours  and  22  minutes 
instead  of  4  hours  and  42  min  after  the  earthquake),  and  overestimates  the  amplitude 
of  the  first  wave  (85  instead  of  70  centimeters).  The  earlier  arrival  time  can  be 
explained  again  because  of  the  no-flux  boundary  conditions  used  along  the  coastlines. 
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Latitude 


Tuticorin  -  Pos.:  08.75N  078. 20E 


Figure  40.  Tuticorin  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 


Tide  gauge  record  Tuticorin  vs.  Model 


Figure  41.  Comparison  of  arrival  times  -  Tuticorin  Tide  Gauge  vs.  Model 
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Latitude 


Kochi  -  Pos.:  09.97N  076. 27E 


Figure  42.  Kochi  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Kochi  vs.  Model 


Figure  43.  Comparison  of  arrival  times  -  Kochi  Tide  Gauge  vs.  Model 
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/.  Colombo  Station  (Sri  Lanka) 

This  station  was  located  at  06.93°1V,  079.83°f7  and  was  approximately 
800  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  44).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  2  minutes  intervals,  was  damaged 
by  the  first  tsunami  wave,  and  did  not  operate  for  5  hours  and  40  minutes  [22], 
According  to  eyewitness  accounts,  the  second  wave  was  the  highest. 

Figure  45  depicts  the  agreement  between  the  tide  gauge  record  and  the 
model  that  the  first  wave  was  positive  (wave  crest)  and  the  highest  wave  was  the 
second.  There  is  a  difference  on  the  arrival  time  (2  hours  and  35  min  instead  of  2 
hours  and  50  min  after  the  earthquake)  for  the  first  wave  and  also  a  difference  in  its 
maximum  amplitude  (80  instead  of  200  centimeters). 

g.  Hanimaadhoo  Station  (Maldives) 

This  station  was  located  at  06. 77°  IV,  073.18°f7  and  was  approximately 
1,200  nautical  miles  from  the  west  side  of  the  north  fault  segment  (Figure  46).  The 
GLOSS  tide  gauge  of  this  station  was  sampled  in  2  minutes  intervals. 

Figure  47  depicts  the  agreement  between  the  tide  gauge  record  and  the 
model  on  that  the  first  wave  was  positive  (wave  crest).  There  is  a  small  difference  on 
arrival  time  (3  hours  and  26  minutes  instead  of  3  hours  31  minutes  after  the  earth¬ 
quake)  for  the  first  wave  and  also  a  difference  in  its  maximum  amplitude  (60  instead 
of  160  centimeters). 
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Latitude 


Colombo  -  Pos.:  06.93N  079.83E 


Figure  44.  Colombo  Tide  Gauge  Position 


67 


Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Colombo  vs.  Model 


Figure  45.  Comparison  of  arrival  times  -  Colombo  Tide  Gauge  vs.  Model 


Latitude 


Hanimaadhoo  -  Pos.:  06.77N  073. 18E 


Figure  46.  Hanimaadhoo  Tide  Gauge  Position 
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Sea  surface  elevation  [cm] 


Arrival  Time 

Tide  gauge  record  Hanimaadhoo  vs.  Model 


Figure  47.  Comparison  of  Arrival  Times  Hanimaadhoo  Tide  Gauge  vs.  Model 
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A  summary  of  simulation’s  results  for  the  tide  gauge  stations  located 
at  South  -  Southeast  of  India  is  presented  in  Table  VII. 


No 

Tide  Gauge 
Station/ 

Field  measurement 

Coordinates 

Arrival  Time 
Difference 
(min) 

Max  wave’s  height 
Relative  error 

1 

Paradip 

20.26°A 

086.70°F 

0 

-0.8 

2 

Vishakhapatham 

17.65°A 

083.28°F 

-10 

-0.4 

3 

Chennai 

13.10°A 

080.32°F 

0 

0.78 

4 

Tuticorin 

08.75°V 

078.20°F 

-10 

-0.44 

5 

Kochi  (Cochin) 

09.97°1V 

076.27°F 

-20 

0.38 

6 

Colombo 

06.93°A 

079.83°F 

-15 

-0.6 

7 

Hanimaadhoo 

06.77°A 

-5 

-0.6 

073.18°F 

Table  VII.  Summary  of  Results  for  the  Tide  Gauge  Stations  Located  at  the  South- 
Southeast  of  India  regions 


Comment: 

(1)  The  relative  error  is  evaluated  as  follows:  Error  ~  hmodei-KbS 
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B.  COMPARISON  TO  SATELLITE  ALTIMETRY 
1.  Jason  1 

One  of  the  three  satellites  used  in  this  study  was  Jason— 1.  This  satellite 
crossed  the  equator  at  2:55  a.m.  UTC,  which  was  two  hours  after  the  earthquake 
[26],  and  its  altimeter  recorded  the  differences  in  the  sea  surface  elevation  in  the 
Bay  of  Bengal.  Figure  48  presents  the  model’s  sea  surface  elevation  two  hours  after 
the  earthquake.  The  black  line  in  this  figure  represents  the  points  used  to  compare 
sea  surface  elevation  in  Figure  49  between  the  satellite’s  altimetry  data  (blue  line) 
and  the  model’s  data  (red  dots).  The  model  estimates  very  accurately  the  leading 
wave  crest  at  about  5°S  and  the  double  peak  structure  between  this  and  the  equator, 
although  the  heights  are  lower  than  the  observations.  However,  between  5°N  to  12°N 
the  model’s  crest  is  in  contrast  to  the  observation’s  trough.  Finally,  between  12°N  to 
20°N  the  model  and  satellite  data  are  in  total  agreement  and  specifically,  the  model 
estimates  very  accurately  the  leading  wave  crest  at  about  20°N. 


Sea  surface  elevation  2  hours  after  the  earthquake 


30°E  45°E  60°E  75°E  90°E  105°E 

Longitude 


Figure  48.  Model’s  Sea  Surface  Elevation  Two  Hours  After  the  Earthquake 
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Jason  1  vs.  Model 


[2  hrs  00  min  after  the  earthquake] 


Latitude 


Figure  49.  Comparison  of  Jason  1  vs.  Model 
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2.  Topex/Poseidon 

The  second  satellite  used  in  this  study  was  Topex/Poseidon.  This  satellite 
crossed  the  equator  at  3:00  a.m.  UTC.  Figure  50  presents  the  model’s  sea  surface 
elevation  two  hours  and  five  minutes  after  the  earthquake.  Again,  the  black  line 
represents  the  points  used  for  comparing  sea  surface  elevation  in  Figure  51  between 
the  satellite’s  altimetry  data  (blue  line)  and  the  model’s  data  (red  dots).  In  this 
case,  the  satellite’s  altimetry  data  contain  several  gaps.  The  model  estimates  very 
accurately  the  leading  wave  crest  at  about  5°S.  The  gaps  in  the  satellite’s  data  between 
5°N  and  the  equator  do  not  show  the  double  peak  structure  that  was  evident  in  Jason- 
1  data.  Between  12°N  to  20°N  the  model  and  the  satellite’s  data  are  in  exceptional 
agreement,  and  in  particular,  the  model  estimates  very  accurately  the  leading  wave 
crest  at  about  20°N. 


Sea  surface  elevation  2  hours  05  minutes  after  the  earthquake 
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Figure  50.  Model’s  Sea  Surface  Elevation  Two  Hours  and  Five  Minutes  After  the 
Earthquake 
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Sea  surface  elevation  [cm] 


Topex  /  Poseidon  vs.  Model 
[2  hrs  05  min  after  earthquake] 


Latitude 


Figure  51.  Comparison  of  Topex/Poseidon  vs.  Model 
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3.  Envisat 

Envisat’s  altimetry  data  were  also  used  in  this  study.  Envisat  crossed  the 
equator  at  4:14  a.m.  UTC.  Figure  52  presents  the  model’s  sea  surface  elevation 
three  hours  and  fifteen  minutes  after  the  earthquake.  The  black  line,  as  previously, 
represents  the  points  used  to  compare  the  sea  surface  elevation  compared  in  Figure 
(53)  between  the  satellite’s  altimetry  data  (blue  line)  and  the  model’s  data  (red 
dots).  In  general,  this  comparison  indicates  that  the  simulation  is  very  accurate  in 
representing  tsunami  waves.  In  addition  to  the  previous  two  cases,  the  reflected  wave 
from  Sri  Lanka  between  10°N  to  15°N  is  reproduced  also  accurately.  The  height  is 
slightly  lower,  due  to  the  lack  of  inundation  algorithms,  which  affects  the  wave  heights 
by  the  introduction  of  artificially  deep  waters  at  the  shallow  depth  near  coastlines. 


Sea  surface  elevation  3  hours  1 5  minutes  after  the  earthquake 
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Figure  52.  Model’s  Sea  Surface  Elevation  Three  Hours  and  Fifteen  Minutes  After  the 
Earthquake 
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Figure  53.  Comparison  of  Envisat  vs.  Model 
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C.  COMPARISON  TO  FIELD  MEASUREMENTS 
1.  Mercator 

During  the  Indian  Ocean  tsunami,  the  Mercator,  was  anchored  about  1.6  kilo¬ 
meters  off  the  Phuket  coast  (Thailand)  and  specifically  at  the  location  07. 715° N, 
098.28 °E.  Based  on  the  yacht’s  depth  gauge  measurements,  the  first  wave  (trough) 
arrived  at  its  location  at  02:38  a.m.  UTC,  one  hour  and  thirty-nine  minutes  after  the 
earthquake,  while  the  three  main  tsunami  waves  had  trongh-to-crest  wave  heights  of 
6.6,  2.2,  and  5.5  meters.  Figure  54  depicts  the  simulated  elevations  and  also  indicates 
the  observed  arrival  time.  The  model  appears  to  simulate  very  accurately  the  arrival 
time  of  the  first  wave  at  the  Mercator’s  location  and  that  this  was  negative  (wave 
trough),  indicating  the  subsidence  on  the  eastern  side  of  the  initialization  area.  How¬ 
ever,  the  amplitude  is  underestimated.  The  reason  for  that,  by  taking  into  account 
similar  difficulties  of  other  models  [31],  is  most  likely  caused  by  misrepresentation  of 
the  water  depth.  The  simulated  water  depth  (grid  depth)  at  the  Mercator’s  location 
was  45  meters  in  comparison  to  an  observed  of  12  meters. 
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Sea  surface  elevation  [cm] 


Arrival  Time 

"Mercator"s  fishfinder  vs.  Model 


Figure  54.  Comparison  of  arrival  times  -  Mercator’s  fishfinder  vs.  Model 
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D.  SYNOPSIS 

From  the  examination  of  the  foregoing  results  it  is  evident  that: 

1.  The  model  accurately  predicts  the  positive  (wave  crest)  first  wave  at  all  tide 
gauge  stations  positions  due  to  the  uplift  of  the  western  side  of  the  initialization  area, 
and  the  negative  (wave  trough)  at  the  location  of  the  Mercator  due  to  the  subsidence 
on  the  eastern  side. 

2.  The  model  simulates  very  accurately  the  arrival  times  and  the  amplitudes 
of  the  first  tsunami  waves  for  the  tropical  Indian  Ocean  (Male,  Gan,  Diego  Garcia, 
Port  Louis,  Salalah,  and  Pointe  La  Rue)  and  the  southwest  Indian  tide  gauge  stations 
(Okha  and  Mormugao).  All  these  stations  received  tsunami  waves  directly  from  the 
source  area. 

3.  The  differences  in  arrival  times  or  in  maximum  amplitudes  between  simu¬ 
lation  results  and  records  from  the  Indian  tide  gauge  stations  (Paradip,  Vishakhap- 
atham,  Chennai,  Tuticorin,  Kochi,  Colombo  and  Hanimaadhoo  (Maldives))  are  due 
to  the  absence  of  inundation  algorithms.  According  to  [22],  the  tsunami  waves  that 
arrived  at  these  stations,  combined  of  waves  traveling  parallel  to  the  earthquake 
rupture  and  waves  reflected  from  the  southwestern  coasts  of  the  Indonesian  Islands. 
Additionally,  in  the  Hanimaadhoo,  Colombo,  Tuticorin,  and  Kochi  stations’  cases,  the 
tsunami  waves  had  to  pass  over  areas  with  depths  less  than  11  meters.  Due  to  the 
boundary  conditions  used  at  these  simulations  (all  positions  with  depths  less  than  11 
meters  were  treated  as  having  depths  of  11  meters)  the  model  predicts  earlier  arrival 
times. 

4.  The  model  underestimates  the  amplitude  of  the  arrived  tsunami  waves  at 
the  Mercator’s  location  due  to  misrepresentation  of  the  water  depth  by  the  used  mesh 
grid. 

5.  By  taking  into  account  that  the  tsunami  wave  amplitude  squared  is  propor¬ 
tional  to  the  potential  energy,  the  results  depict  the  west  and  southwest  propagation 
of  most  of  the  energy  due  to  the  orientation  of  the  earthquake  rupture  (see  Figure 
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55).  Only  a  fraction  of  the  energy  propagated  toward  the  south  and  southeast,  which 
is  in  agreement  with  [22], 


Maximum  heights  of  waves  in  meters 
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Figure  55.  Energy  Propagation  of  the  Simulated  Tsunami  Waves 


6.  According  to  [22]  the  mid-ocean  ridges  played  a  major  role  as  wave  guides 
that  transferred  the  tsunami  energy  to  far-held  regions  outside  the  source  area  in  the 
Indian  Ocean.  The  model  reproduced  this  feature  very  accurately  as  illustrated  in 
Figure  55  in  combination  with  Figure  1.  Finally,  Figure  55  presents,  the  areas  where 
the  tsunami  focused  most  of  its  energy. 
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V.  CONCLUSIONS  AND 
RECOMMENDATIONS 

This  study,  as  part  of  a  research  effort  regarding  the  development  of  a  triangu¬ 
lar  discontinuous  Galerkin  oceanic  shallow  water  model,  formatted  real  bathymetry 
data  of  the  Indian  Ocean  and  simulated  the  propagation  stage  of  the  Indian  Ocean 
tsunami  that  occurred  on  December  26,  2004.  The  model’s  main  advantage  is  related 
to  the  geometric  (grid)  flexibility,  which  made  it  possible  to  represent  the  tsunami’s 
propagation  with  different  resolutions  throughout  the  Indian  Ocean. 

The  model’s  results  were  compared  to  the  records  of  fifteen  tide  gauge  stations 
around  the  Indian  Ocean,  satellite  altimetry  (Jason-1,  Topex/Poseidon  and  Envisat 
data),  and  field  measurements.  These  results  showed  that  the  model  is  capable  of 
producing  accurate  estimates  of  arrival  times  at  distant  locations.  Some  considerable 
differences  were  noticed  from  the  above  comparisons  due  to  influences  by  reflections, 
misrepresentation  of  the  water  depth  by  the  mesh  grid,  and  the  no-flux  boundary 
conditions  serving  as  a  proxy  for  inundation  algorithms. 

The  thesis  also  derived  the  analytic  solution  of  the  linear  Munk  problem  (two- 
dimensional  problem)  for  future  use  in  validating  the  diffusion  operators  that  were 
not  used  in  the  tsunami  simulations. 

Future  studies  should  involve  the  following: 

1)  The  use  of  inundation  algorithms  in  order  to  complete  the  simulation  of 
the  Indian  Ocean  tsunami  December  26,  2004  (propagation  and  run-up  stages),  by 
using  this  DG  model.  These  results  could  be  compared  against  the  already  existing 
de-tidal  records. 

2)  The  simulation  of  the  Indian  Ocean  tsunami  December  26,  2004  (prop¬ 
agation  and  run-up  stages),  by  using  this  DG  model  and  more  complicated  initial 
conditions  (for  example  more  than  two  planes). 

3)  The  simulation  of  earthquakes  with  different  magnitude  for  the  same  region 
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of  the  Indian  Ocean,  or  the  simulation  in  other  parts  of  the  world,  in  order  to  identify 
the  regions  that  are  likely  to  be  severely  affected  by  tsunami  waves,  and  to  establish 
seismic  criteria  for  issuing  tsunami  warnings  in  the  case  of  actual  tsunamis. 

4)  The  effects  of  storm  surge  modeling  (propagation  and  run-up)  caused  by 
hurricanes,  cyclones  and  typhoons  on  coastal  areas. 
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APPENDIX  A.  DERIVATION  OF  ANALYTIC 
SOLUTION  FOR  THE  MUNK  PROBLEM 


In  Stommel  ’s  theory,  the  inclusion  of  the  bottom  friction  raises  the  order  of 
the  governing  equation  and  enables  the  beta  effect  to  produce  a  western  boundary 
current  [19],  which  makes  Stommel’s  solution  not  quite  realistic.  The  ocean  currents 
are  concentrated  in  the  upper  kilometer  of  the  ocean,  they  are  not  barotropic,  and 
they  are  not  independent  of  depth. 

Walter  Heinrich  Munk  noted  that  the  frictional  effect  due  to  small-scale  eddies 


is  crucially  important  in  ocean  dynamics,  and  solved  the  problem  by  considering  the 
lateral  eddy  viscosity  to  explain  the  ocean  current  distribution,  based  on  the  observed 
wind  system. 

Consider  a  flat-bottomed  ocean  of  depth  H  driven  by  wind  stress  at  its  surface. 
The  momentum  shallow  water  equations  in  two-dimensions,  where  bottom  drag  has 


been  neglected  while  the  lateral  eddy  viscosity  is  included,  are: 
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Pressure  can  be  eliminated  by  cross  differentiating  these  equations  and  taking  the 
difference  J^(Al)  —  J^(A 2)  : 


d  <9 

&(/u)  +  SyW 
du  dv  df 

!{Tx  +  d^,)  +  lh,v 


dv  du.  Id, 

=  "v"(&-^)  +  yadOTr'T) 

,dv  du.  ^d, 

=  vV^dx  “  ft/) +  y  ai<cm'  t) 


(A.3) 


The  use  of  the  continuity  equation  §y  +  =  0  and  beta-plane  approximation 

(3  —  into  Equation  (A.3),  gives  the  following: 
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( curlf ) 
( curlr ) 


(A.4) 


The  Equation  (A. 4),  after  multiplying  by  pQ,  integrating  from  the  bottom  (-H)  to  the 
surface  (zero)  in  the  vertical  with  boundary  conditions  w  —  0  at  z  —  0  and  z  =  —H, 
assuming  zero  bottom  friction,  and  using  the  definition  of  the  streamfuntion: 


gives  the  following: 


Po 


/  udz  —  —  TT 

I-h  ay 

f°  ,  „ ,  <9T 

/  vdz  =  Mv  = 

l—H  OX 


^  = 


Coriolis 


+  curlf0 

Lateral  Friction  Winds 


(A.5) 

(A.6) 

(A. 7) 


According  to  Munk,  in  the  western  boundary  layer  (WBL)  Coriolis  force  is 
being  balanced  by  the  lateral  friction  while  at  the  interior,  Coriolis  is  being  balanced 
by  the  wind’s  effect: 
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Consider  the  idealized  model  for  the  North  Atlantic  subtropical  gyre,  where: 

-  Lx  —  5000 km 

-  Ly  =  2500 km 

-  Lm  is  the  width  of  the  WBL,  (Lm  -C  Lx) 

-  Tox  —  —A  cos  j1  is  the  x-component  of  the  wind  force 

-  Toy  =  0  is  the  y-component  of  the  wind  force 

-  A  =  0.1^. 
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i.  Interior  Solution 

The  solution  of  the  equation 
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is  as  follows: 
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Thus,  at  the  interior,  Lm  <  x  <  Lx  and  0  <  y  <  Ly,  the  streamfunction  solution  can 
be  written  as 


V(x,y)  =  (Lx  -a:)  sin  (^)  (A.8) 

The  use  of  Equation  (A.8)  and  Equations  (A. 5)  and  (A. 6)  of  the  streamfunction’s 
definition  gives 
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The  momentum  shallow  water  equations  in  two-dimensions  can  also  be  written  as 
follows: 
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where  h(x,  y)  is  the  free  surface  elevation.  The  substitution  of  the  Equations  (A. 9), (A.  10) 
into  (A.  11), (A.  12)  gives  the  following: 


A  fir  (iry\  Auir2  ( irys 

hMix' v)  =  (L‘  -x)sm{Tj-  h  cos  U, 


+  C1!  (A.  13) 


where  C\  is  an  unknown  constant,  and  multiplying  by  the  gravitational  acceleration 
g,  gives  the  geopotential  height  solution: 

Auir2 


Mx'v)  =  (L*-X)sm  -  MPJ. I 


cos 


7 ry 


A  gC i  (A. 14) 
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Ai  =  0 


Since  the  WBL  solution  cannot  grow  exponentially,  it  must  be  that  C2  =  0,  and 
the  three  remaining  unknowns  Ci,  C3,  and  C4  can  be  evaluated  by  using  the  following 
conditions: 


(1)  ^U=o  =  0  =►  ^|x=0  =  0 

(2)  v\x=q  =  0  fU=0  =  0  =>  |||a;=0  =  0 

(3)  at  distance  Lm  the  two  solutions  (WBL  and  interior)  must  be  the  same 
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The  use  of  the  above  three  conditions  gives  the  following: 


T  ( x )  = 


AnLx 

BL 


y  L 


x  /  VSx  1  VSx ' 

1  —  e  2L ™  cos  — - 1 - 7=  sin  — - — 

V  2  Lm  VS  2  Lm/ 


Thus,  at  the  WBL,  0  <  x  <  Lm  and  0  <  y  <  Ly,  the  streamfunction  solution  can  be 
written  as 


T  (x,y)  = 


AnL a 


X  (  VSx  1  VSx ' 

1  —  e  2Lm  [  cos  77-7 - 1 - 7=  sm 
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VS  2  Lr 


sin  (m  (A. 15) 


As  for  the  interior  solution,  the  use  of  Equation  (A. 15)  and  Equations  (A. 5)  and  (A. 6) 
of  the  streamfunction’s  definition  gives  the  following: 

Att2Lx  f  _ sl_  (  VSx  VS  .  VSx\\  (ny\ 
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3  p0H  (3LyLm  2  Lm  \Ly ) 

With  the  same  way  as  for  the  interior  solution,  the  substitution  of  the  Equations 
(A. 16), (A. 17)  into  (A. 11), (A. 12)  gives 
AnfLx 
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where  C2  is  an  unknown  constant,  and  multiplying  by  the  gravitational  acceleration 
g  gives  the  geopotential  height  solution: 
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